import math
import numpy as np
import matplotlib.pyplot as plt

# 西瓜数据集3.0alpha
# 密度 density
d = np.array([0.697, 0.774, 0.634, 0.608, 0.556, 0.403, 0.481, 0.437,
              0.666, 0.243, 0.245, 0.343, 0.639, 0.657, 0.360, 0.593, 0.719])
# 含糖率 sugar content
s = np.array([0.460, 0.376, 0.264, 0.318, 0.215, 0.237, 0.149, 0.211,
              0.091, 0.267, 0.057, 0.099, 0.161, 0.198, 0.370, 0.042, 0.103])
# 1 表示好瓜 0 表示坏瓜
label = np.array([1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0])

# 变量
step = 0 # 迭代步数
old_l = 0 # 上一次迭代中似然函数的值
beta = np.zeros(3) # beta = (w;b) 权值和偏置

# 迭代循环（牛顿法）
while 1:
    cur_l = 0 # 当前对数似然函数（公式3.27）
    beta_x = np.zeros(17) # 17个样本的beta^{T}\hat{x} = w^{T}x + b

    # 计算似然函数l(beta)（公式3.27）
    for i in range(17):
        beta_x[i] = (np.dot(beta, np.array([[d[i]], [s[i]], [1]])))[0]
        cur_l = cur_l + ((-label[i] * beta_x[i]) + math.log(1 + math.exp(beta_x[i])))

    if abs(cur_l - old_l) < 0.001:
        # 如果似然函数的值无变化，则说明迭代收敛，停止迭代
        break

    # 否则进行一次迭代，更新beta的值
    step += 1
    old_l = cur_l
    d1 = 0 # l关于beta的一阶导数
    d2 = 0 # l关于beta的二阶导数
    p1 = np.zeros(17) # 17个样本在公式3.30和3.31中p1(\hat{x}_{i};beta)的值

    # 计算一阶导数和二阶导数
    for i in range(17):
        p1[i] = 1 - 1 / (1 + math.exp(beta_x[i]))
        x_i = np.array([[d[i]], [s[i]], [1]])
        d1 = d1 - np.dot(label[i] - p1[i], x_i) # 公式3.30
        d2 = d2 + p1[i] * (1 - p1[i]) * np.dot(x_i, x_i.T) # 公式3.31

    # 更新权重beta
    beta = beta - (np.dot(np.linalg.inv(d2), d1)).T # 公式3.29

# 测试
predict = [] # 对率回归的预测结果
acc = 0 # 预测正确的数量
for i in range(17):
    logit = (np.dot(beta, np.array([[d[i]], [s[i]], [1]])))[0] # 公式3.19、3.21

    if (logit >= 0 and label[i] == 1) or (logit < 0 and label[i] == 0):
        # 预测正确
        acc += 1
    
    predict.append(1 if logit >= 0 else 0) # 保存模型预测的结果

print('Accuracy : ', float(acc) / float(len(label))) # 输出正确率
print('Label    : ', list(label)) # 输出标签
print('Predict  : ', predict) # 输出预测结果

# 绘图
x = np.arange(0.1, 0.9, 0.1)
# 分界线为：0 = beta0 * d + beta1 * s + beta2 * 1
# 移项：s = (-beta0 * d - beta2) / beta1
y = np.array((-beta[0][0] * x - beta[0][2]) / beta[0][1])
good = plt.scatter(d[0:8], s[0:8], color = 'g', marker = 'o', label = 'good')
bad = plt.scatter(d[8:17], s[8:17], color = 'r', marker = 'x', label = 'bad')
plt.plot(x, y)
plt.xlabel('density')
plt.ylabel('sugar content')
plt.legend()
plt.savefig('./3_3.png')
